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Abstract 

We present an NLO analysis of the massive vector current correlator at temperatures above 
a few hundred MeV. The physics of this correlator originates from a transport peak, related 
to heavy quark diffusion, and from the quark-antiquark threshold, related to quarkonium 
physics. In the bottom case both can be studied with separate effective theories, but for charm 
these may not be accurate, so a study within the full theory is needed. Working in imaginary 
time, the NLO correlator can be computed in unresummed perturbation theory; comparing 
with lattice data, we find good agreement. Subsequently we inspect how non-perturbative 
modifications of the transport peak would affect the imaginary-time correlator. The massive 
NLO quark-number susceptibility is also contrasted with numerical measurement. 
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1. Introduction 



Heavy (charm and bottom) quarks are excellent probes for the properties of the hot QCD 
plasma generated in heavy ion collision experiments. On the theoretical side, the existence of 
a mass scale M ^ vrT S> 200 MeV renders the heavy quarks relatively tractable, permitting 
for an interpolation between the simple static dynamics of the infinite mass limit and the 
high mobility case manifested by lighter quarks. It is particularly fortunate that two heavy 
flavours are available, offering for a handle on the functional dependence on M/ttT. On the 
experimental side, heavy quarks and quarkonia are readily tagged because of their distinctive 
leptonic decays. Indeed thermal modifications of the bottomonia spectra were among the 
first spectacular results produced by the LHC heavy ion program [1]. 

For the bottom quark case, recent years have seen significant progress in theoretical studies 
of the main phenomena involved, namely single quark "transport" (diffusion, kinetic equi- 
libration) as well as physics near the quark-antiquark threshold (quarkonium dissociation, 
chemical equilibration) (cf. rcf. [2] for a review and refs. [3]-[12] for recent contributions). 
Largely this progress has been achieved through the use of modern effective field theory meth- 
ods (Heavy Quark Effective Theory, or HQET, for single quarks; Non-Relativistic QCD, or 
NRQCD, for quarkonium physics). Once properly formulated the effective field theory ob- 
servables can be measured non-perturbatively with lattice Monte Carlo methods, and indeed 
first results suggest that these avenues may lead to substantial progress [13]-[17]. 

In the charm quark case, however, it is not guaranteed that the heavy quark expansion 
converges fast enough to yield quantitatively accurate results for all observablcs of interest. 
On the other hand, it no longer appears prohibitively expensive to treat charm quarks as 
"light" degrees of freedom in lattice QCD. Indeed, results have appeared concerning both 
thermodynamic quantities [18, 19, 20] and imaginary-time correlators relevant for determin- 
ing dynamical properties of the system [21] (earlier works can be found in refs. [22]-[25] and 
references therein). Yet, a modest scale hierarchy between ttT and M does exist, and con- 
sequently systematic errors, both from lattice artifacts and from the unavoidable analytic 
continuation [26], are likely to be harder to control than in the light quark case. In fact, 
given that systematic errors related to analytic continuation remain substantial even for light 
quarks [27], further crosschecks appear welcome. 

The goal of the current study is to derive results within next-to-leading order (NLO) 
perturbation theory which may help in the analysis of lattice data such as those in ref. [21]. 
In order to allow for a direct appreciation of Euclidean measurements, we inspect how specific 
modifications of transport properties and threshold features manifest themselves in imaginary 
time. Ultimately, in accordance with the philosophy of ref. [28] and practical tests of refs. [27, 
29], the goal would be to subtract "trivial" ultraviolet features from continuum-extrapolated 
lattice data, in order to allow for a model- independent extraction of real-time physics [30]. 
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On the technical side, the current work represents a continuation of our earher study [31], 
in which the massive vector current spectral function was computed at NLO in the domain 
M S> ttT (keeping only those thermal effects which are not exponentially suppressed). Here 
we keep the full mass dependence, permitting for an extrapolation also to the regime M <C vrT, 
as well as contributions from the transport peak at u <^ 2M which were omitted in ref. [31]. 
Moreover we work directly in imaginary time, which has the benefit that the usual problems 
of convergence at small ui are milder. In particular, at NLO infrared safe results can be 
obtained without resummations, similarly to what has previously been achieved for gluonic 
observables [32, 29]. 

The plan of this paper is the following. After specifying the observables considered and 

discussing the methods employed (sec. 2), we outline the qualitative structure of our find- 
ings in sec. 3. The detailed analytic and numerical results of the strict NLO analysis com- 
prise sec. 4, whereas in sees. 5 and 6 the effects of non-perturbative modifications of the 
transport peak and quarkonium threshold, respectively, are inspected. Sec. 7 presents our 
conclusions; appendix A results for all the "master" sum-integrals considered; and appendix B 
details related to renormalization. 



Observables and methods 



2.1. Basic definitions 



The main quantity considered is the vector current correlator related to a massive flavour. 
Like in lattice QCD we work in Euclidean signature, with the usual thermal boundary con- 
ditions imposed across the time direction. Then the correlator is defined as 

d 

(2.1) 

//=0' 



d „ 

G^{t) ^ ^ /((Vi7MV')(r,x)(Vi7^V)(0,0) 



= C!qq{t) - G.i.i{t) 



0<T< ^ , 



(2.2) 



where the Dirac matrices arc Minkowskian, and a sum over spatial indices is implied. Be- 
cause of current conservation the charge correlator is independent of r, and we denote the 
corresponding "susceptibility" by 

X ^ iyrG,,{T) = pGooiO) , /3 ^ i . (2.3) 

For future reference let us also record the free massless results for these correlators [33] : 
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Here Ca = Nc = 3 refers to the number of colours. Later on the group theory factor 
Cp = {N^ — l)/{2Nc) win also appear. Spacetime dimension is denoted hy D = d+1 = 4 — 2e. 
At NLO, we find it convenient to compute the correlators and Gqq. The spatial part 
is then obtained from eq. (2.2). Verifying explicitly the r-independence of Gqo provides 
for a nice crosscheck of the computation. 



(2.6) 



2.2. Wick contractions for the vector current correlator 

The Wick contractions for are as given in ref. [31]. Denoting 

Q = (a;„,0), Ap = P^ + M\ 
the leading-order (LO) vector correlator reads, in momentum space, 

Here {P} denotes fermionic Matsubara momenta, and f^^py = TY2{p„} Jp- 

At NLO, we have to decide on a meaning of the renormalized mass. Although conceptually 
subtle, it is technically convenient to employ a pole mass; then the bare mass parameter, M|, 
can be expressed as M| = + where at NLO 
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M2 ^ 3 



Here ii is the scale parameter of the MS scheme, terms of 0{e) were omitted, and 
efc = |k| , Ep = ^p^ + M2 , Epk = v/(p-k)2 + M2 . 



(2.10) 



(2.11) 



Because of Lorentz invariance the vector P in eq. (2.8) can be chosen at will, as long as we 
set pq = ziziEp after carrying out the X-integral; this means that the vector p in eq. (2.9) is 
arbitrary. The corresponding counterterm contribution reads 



4Ca(5M2 



f 

{P} 



1^ 



+ 



4M2 -{D- 2)Q2 

^l^P-Q 



(2.12) 



and it is often convenient to identify p of eq. (2.9) as the integration variable of eq. (2.12). 
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The "genuine" 2-loop graphs amount to 



■^K{P} ^ 

(D-2)2 {D-2f 2(1) -2) 4(L»-2)M2 
2{D - 2) A{D - 2)M2 - {D- 2)2q2 2{D - 



K2ApAp_Q ^'A2,Ap_Q i^2A^_^Ap_Q 

_^ 4(1) - 2) 4(L> - 2)M2 - (L> - 2)2q2 

ApAp_j^Ap_Q ApAp_^.Ap_g 

16M2 + 2(L» - 2)2iv: • Q - 4(Z) - 2)Q2 igM"^ - 4(Z) - 2)Q^M^ 
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2.3. Wick contractions for the susceptibility 

In the case of the zero components, viz. Gqq, the LO correlator reads, in momentum space, 
whereas the counterterm graph can be expressed as 



^C^SM^t i^-^^ (2-15) 

{■P} P P-Q P 



The genuine 2-loop graphs amount to^ 
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D-2 D- 2 2 4M2 



+ ^^. + 
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^The appearance of the "energy variables" in the numerators implies that there is a certain redundancy in 
the basis: 

[ I Q'^l ] (216) 

T^lP^X K^^pAp_j,Ap_Q ApAp_^.Ap_QAp_^_Q X^ApAp.^Ap.QAp.^,.^ / ' 
We have used this to eliminate terms containing Q^ei in the numerator. 



2{D-2) ^ (£'-2)(4E2 + Q2) 
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^ 4(17-2) iD-2)i-iEf,+Q^) 
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2.4. Matsubara sums and spatial integrals 

The next step is to convert the momentum-space expressions to configuration space. If we 
denote the above result (eqs. (2.7), (2.12), (2.13)) by (5y(a;„), then the conversion is obtained 
as 

G^(r)=T^e-^-"^G^(a;„), (2.18) 

and similarly for Gqq. At NLO we are thereby faced with a three-fold Matsubara sum. Making 
use of standard techniques, reviewed in some detail in ref. [31], these sums can be carried 
out in a closed form, whereby we are left with integrals over at most two spatial momenta 
(i.e. two radial directions and one angle). Intermediate results at this stage are displayed for 
all individual master sum-integrals in appendix A, and for their sums in appendix B, in the 
latter case with the cancellation of 1/e-divergences verified as well. 

Two interesting crosschecks are available. First of all, all r-dependent terms disappear 
from Gqq, cf. eq. (B.3). Second, individual parts of the expressions contain "contact terms" 
oc ^^(r), where (5^ denotes the periodic Dirac-delta. These arise from structures that are 
independent of a;„, but also from sum-integrals in which appears in the numerator. It can 
be verified, however, that all contact terms cancel, both at LO and at NLO. 

It remains to carry out the spatial integrals. We write 



/ -I /• 

v/p,k J p,k J z 



(2.19) 



where the normalization of the angular variable z = cos^p^k is chosen so that 1 = 1. 
Depending on the case, it is convenient to substitute variables as 

r 1 rKk 1 r^tk 

X ^ ^ Je- ^ ^ Je- "^'^''^^ ' ^^"^^^ 
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where = ± fc)^ + M^, e^^ = |p± /c|. The angular integrals are doable in most cases. 
In addition, it is also possible to carry out partial integrations with respect to the radial 
directions, which helps to reduce the number of independent terms (for the massless case, see 
ref . [34] for a recent discussion) . Closed massless loop integrals are typically solvable, but in 
general the integrations remain to be carried out numerically. Our final expressions are given 
in the next two sections, cf. eqs. (3.2)-(3.5), (4.1), (4.4), (4.5). 



3. Leading order results and qualitative pattern 

In order to illustrate the qualitative structure of the results, we recall in this section the LO 
expressions for the quantities considered. In general, two types of contributions appear: those 
that depend on r, and those that are constant. To display the r-dependence we introduce 
the periodic dimensionless function 

(Bi+---+Bfc)(^-r)+(Bfc+i+-+B„)r I (Bi+-+Efc)r+(Bfc+i+-+B„)(/3-T) 

^E,-E, {^)- [e^^^i ± 1] . . . ± 1] • ^""^^ 

In the denominator, the sign is chosen according to whether the particle is a boson or a fermion 
(at NLO, there is only one boson, with the on-shell energy denoted by ek, cf. eq. (2.11)). We 
also denote £'2Ep = ^EpEp- With these conventions, the LO result for Gy reads 

G^°WI.-aep. = -Ginr)U.. = -'^C^f{^ + ^)D,Epir), (3.2) 



Gt°WI„.. - -4C. / . ' (3.3) 



i 



Here rip denotes the Fermi distribution (ne denotes the Bose distribution) . The susceptibility 
only contains a r-independent part: 

Tx"° = GJi? = -4Ca / TnUEp) . (3.4) 
Jp 

Finally, the constant part of the spatial correlator is obtained through the use of eq. (2.2): 

'p ^p 



Gr(r)|_, = -^Ca / (l - -^)TnUE,) . (3.5) 



To our knowledge none of these leading-order expressions can be integrated in terms of 
standard elementary functions. 

In terms of the spectral function, viz. Pjj(w), the constant contribution in eq. (3.5) arises 
from (an infinitely narrow) transport peak around u = 0, whereas the "fast" r-dependence in 
eq. (3.2) originates from the quark-antiquark continuum at > 2M. Around the middle of 
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the Euclidean time interval, D2^^{(3/2) = —2Tn'^{Ep). Therefore, both terms contribute in 
a comparable manner at large Euclidean time separations, even though they originate from 
completely different types of physics (see also ref. [35]). 

At NLO, the same three structures appear as at LO: a constant Gqq, as well as a spatial 
correlator G^^ which has both a constant and a r-dependent part. It is generally believed that 
the perturbativc series for the constant part of G^^ breaks down at some point and that in the 
full theory, Gjj has no projection to the Matsubara zero mode; the reason is that the spatial 
components of the vector current do not couple to conserved charges. The corresponding 
"slow" T-dependence then reflects the physics of a "smeared" transport peak. Yet it remains 
true that single quark physics from the transport peak and quark-antiquark physics from 
the threshold region are expected to contribute in a comparable manner to the Euclidean 
correlator around r = ^/2. (This is demonstrated explicitly in figs. 3, 4 below.) 



NLO results 



4.1. Susceptibility 

Proceeding to NLO, the result for the susceptibility is obtained from eq. (B.3) after partial 
integrations: 



Jp Jk . 



. E,. 



1 - 



■ (4.1) 



This can be shown to agree with what can be extracted from ref. [36] as Tc?^ p(T, /Lt)|^^Q. 
(The substance of the information is already there in ref. [37].) The massless loop evaluates 
to Jj, nji{ek)/ef, = r^/12, and in the limit M S> ttT its effect can be interpreted as an effective 
mass correction to the LO result of eq. (3.4), M^+g^T^Cp/G [38]. On the other hand, 

in the massless limit M <^ ttT eq. (4.1) reduces to the well-known correction (cf. [39] and 
references therein) 

A numerical evaluation of eqs. (3.4), (4.1) is shown in fig. 1. For the gauge coupling we have 
inserted the 2- loop value for the thermal coupling g^/T from ref. [40], and for comparison with 
quenched lattice data from refs. [41, 21] we assume that Tc/Aj^ 2± 1.25, T/Tc c± 1.45. In the 
massive case, where no continuum extrapolation is available in ref. [21], we choose to compare 
with results obtained on a lattice 128^ x 48. The quark mass cited, mMs(mc) = 1.094(1) GeV 
is phenomenologically on the low side [42]; based on eq. (2.10), viz. 



M ~ mMs{mc) { 1 + 



(47r)2 



+ Oig' 



(4.3) 
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Figure 1: The quark-number susceptibility, eqs. (3.4), (4.1), normalized to the free result from 
eq. (2.5). The lattice result at M = comes from ref. [41], that at M > from ref. [21]; both are 
quenched (iVf = 0) and only the former represents the continuum limit. The uncertainty of M/T is 
our estimate (cf. the text). Lattice results for charm by other groups can be found in refs. [19, 20]. 



we get M ~ 1.3 GeV but assign a large error to this number. Ref. [21] also cites Tcjyfa = 
0.630(5), ^/a = 428 MeV, so we estimate M/T ~ 3.3 it 0.5, but with substantial systematic 
uncertainties, from lattice artifacts, string tension measurement [43], quenching, as well as 
perturbative input. With this in mind, the excellent agreement seen in fig. 1 is remarkable, 
and supports the long-held belief that all quarks, and heavy quarks in particular, are well 
described by the weak-coupling expansion at surprisingly low temperatures. (In the massless 
case the good agreement is consistent with the recent unquenched study of ref. [44] , in which 
a similar "dimensionally reduced" resummation scheme was used as here [36], however an 
almost perfect match already at NLO may be somewhat coincidental.) 



4.2. Vector current correlator 

The vector current correlator is obtained from eqs. (B.l), (B.2) after partial integrations. Its 
r-dependent part reads 
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(4.4) 



where P denotes a principal value, and A^^ 
contribution reads 



etc are defined in eq. (A. 2). The constant 
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(4.5) 



The spatial correlator subsequently follows from eq. (2.2), with the susceptibility inserted 
from eq. (4.1). 

Although the expression in eq. (4.4) is finite, its numerical evaluation is non-trivial. There 
are at least three separate challenges: at small A; various parts of the expression are divergent, 
and care must be taken in order to avoid significance loss in their cancellation; at A; = p there 
is a pole which is defined in the sense of a principal value; and at large A; there is a vacuum 
part which decreases only slowly (although it is integrable). For the reader's benefit, let us 
briefly specify how we have dealt with these challenges. 

• The small-A; divergence originates from the terms integrated over z in eq. (4.4) and 
from terms where the integral had already been carried out. For the latter type the 
divergent part reads 
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containing both vacuum and Bose-enhanced structures. The integral needs to be 
carried out precisely enough such that the cancellation takes duly place. 

• The principal value integration can be handled for instance by reflecting the range p > k 
into the range p < k: 

f°° /"^ r 1 /k\ 1 

/ dpp'^^{p) = k^ dx x'^(l){kx) + ^(j)(-) . (4.7) 
Jo Jo X ^ 

Here has to be evaluated precisely enough for cancellations at a; = 1 to take place. 

• A possible way to accelerate the convergence at large k is with the help of the function 
6{k) in eq. (4.4). (Note that a power tail only appears in the vacuum part.) The 
simplest subtraction removes just the leading asymptotic behaviour —SEp/k"^, e.g. 

but of course more refined choices can be envisaged. (By replacing A by a "gluon 
mass" it would be possible to take fcmin — > and still carry out J^dk9{k)/k exactly, 
cf. eqs. (B.22), (B.23), but the price to pay is that then the small-A; range has more 
structure than before, with a would-be divergence only cut off at < A.) 

For a transparent representation of the numerical results, we consider two different nor- 
malizations. A simple and theoretically clean possibility is to normalize the results to the 
free massless expression from eq. (2.4). Another reference point is to make use of the full 
NLO spectral function in vacuum [45]-[48] : 

P^H = -n. - 2M) ^-'-' - '^'^^ + '"'^ + - 2M)ggg# ( 



+uj{uj^ - 4M2 



9 w(a;^-4M2) 3, o 



{u' + In ^^^^3 - g (c^^ + 6M^) 



(4.9) 

where the function L2 is defined as 

L2{x) = 4Li2(x) + 2Li2(-x) + [21n(l-x) + ln(l + a:)]lnx. (4.10) 

There is no transport peak in the vacuum expression, and recalling eq. (2.2), the correspond- 
ing spatial correlator can be obtained through 

poo J cosh ( f — r ) a; 

GTir) = / —{-pT){^) } . (4.11) 

Jo smh ^ 
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xT xT 

Figure 2: Left: The vector correlator, normalized to the free result from eq. (2.4), for M/T = 
1.0,2.0,3.0,4.0,5.0 (top to bottom). Right: The vector correlator, normalized to the "reconstructed" 
result from eqs. (4.9), (4.11), for the same masses (left to right). The lattice results are from ref. [21]; 
they are quenched (A^f — 0) and do not contain a continuum extrapolation [N^ — 48, Ns — 128). 

A normalization with respect to a similar "reconstructed" correlator G^^''^{t) has been used 
in ref. [21] (see also ref. [49]), and may be useful for phenomenological purposes, although 
from the theoretical perspective it induces new systematic uncertainties. 

In fig. 2 we show our results in both normalizations, compared with lattice data from 
ref. [21]. (The gauge coupling has been fixed as explained in connection with fig. 1; at 
very small r a different running would appear reasonable but in the absence of an NNLO 
computation and continuum-extrapolated lattice data, we stick to the simplest choice in the 
following.) Like in fig. 1, an excellent agreement is found at large time separations, if a value 
M/T ~ 3.5 is assumed. (The discrepancy in fig. 2(left) at small r is probably due to the 
missing continuum extrapolation.) 



5. Modification of the transport peak 

As mentioned in sec. 3, the correlator G^^ has a constant (r-independent) part at LO and at 
NLO, but within the full dynamics this is expected to turn into a slowly evolving function. 
The purpose of this section is to estimate how precisely Euclidean data should be measured 
in order to resolve the slow time dependence. 
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Figure 3: Left: The effect of a modified transport peak, for M/T = 3.5 and 2'kTD = 1.0, 2.0, 3.0, 
4.0, 5.0 (cf. eq. (5.1)). Right: A magnification of the large-r region. It is clear that a high precision is 
needed for resolving the diffusion coefficient from the massive vector current correlator. (The curves 
could be put on top of lattice data through a minor change of M/T, but we have refrained from doing 
this in the absence of a continuum extrapolation of the reconstructed correlator.) 



In order to reach this goal, we model the transport peak through a Lorentzian shape, 

pfXu:) = 3Dx . — rriT^ ■ (5-1) 
^ ^ oj^ + r/2 cosh(2:^) ^ ^ 

Here D corresponds to the heavy flavour diffusion coefficient. The Lorentzian shape can be 
correct only at small frequencies, |a;| <^ vrT, cf. e.g. ref. [50]; we have chosen to cut it off at 
large frequencies through the same recipe that has been used in the massless case [27]. The 
susceptibility x is fixed according to ref. [21], x/T'^ = 0.20894(1). We then vary D, in each 
case tuning the other parameter rj so as to keep the area under the transport peak fixed at 
the value predicted by the NLO expression, eqs. (3.5), (4.1), (4.5), (2.2): 

' ^'^G^ir) = r^^4^ = [G'I;° + G'r]_, . (5.2) 



/3 70 " ' ' Jo IT ^00 

Subsequently the "correct" G,-,-{t) is obtained by replacing the part [G^^ + G^,^°]^^^^^ through 
a r-dependent function, G^^\t), determined by p^^^ via eq. (4.11). As a guideline, we recall 
that 2itTD ~ 3 — 5 may be expected according to refs. [14, 15]; that ref. [21] found even 
2'kTD ~ 2; and that in the massless case values down to 2ttTD ~ 1 appear possible [27]. 

Results are shown in fig. 3. We observe that if 2'kTD ~ 1 — 2, then there is hope of 
resolving it with high enough precision of the lattice data (it appears that statistical errors 
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Figure 4: Left: The effect of a modified amplitude of a resonance peak, for M/T = 3.5, AM/M = 0.1, 
and A = 0.1, 0.2, 0.3, 0.4, 0.5 (cf. discussion below eq. (6.1)). Right: The effect of a modified shape 
of a resonance peak, for A = 0.5 and AM/M = 0.02, 0.04, 0.06, 0.08, 0.10. 



should probably be reduced to 20% of the current ones to be sensitive to the features of the 
transport peak; obviously improvements in statistical accuracy need to be accompanied by 
a corresponding decrease in systematic uncertainties). In the case 2ttTD ~ 3 — 5 a reliable 
determination from the massive vector current correlator appears challenging. 



6. Modification of the threshold region 



The second qualitative structure affecting the massive vector current correlator, the quarko- 
nium threshold region inducing a "fast" time dependence from the energy scale ~ 2M, is 
also expected to undergo drastic changes in the interacting theory. At low temperatures, the 
spectral function is characterized by quarkonium resonances; at high temperatures, these are 
expected to move, broaden, and eventually dissolve into a mere threshold enhancement. 

Based on our earlier investigations [51, 52], we expect that in the temperature range of 
interest there is at most one resonance peak in the vector channel spectral function, placed 
slightly to the left from the free quark-antiquark threshold. Akin to eq. (5.1), we model this 
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by a skewed Breit-Wigner shape, ^ 



^rV).^^. "'^-2M.A«, (6.1) 



constructed so as not to contribute to the transport coefficient. To reduce the number of 
free parameters to two, we set AM = 27 in the following. The contribution from such a 
peak, determined through eq. (4.11), is added to the thermal NLO result as such, and to 
the vacuum result of eq. (4.9) with A — t- 5A, 7 — )■ 7/5, keeping the area under the peak 
roughly invariant. Obviously these choices are arbitrary, but they should nevertheless convey 
a qualitative impression on the importance of resonance contributions. Based on ref. [31], 
in which a resonance peak around a threshold was matched to the thermal NLO spectral 
function above the threshold, we expect that A<0.5 and AM/M<0.1. 

Results are shown in fig. 4. Comparing with fig. 3, it can be seen that a change of the 
spectral shape around the threshold region affects the Euclidean correlator at smaller r than 
a change of the transport peak. In fact, if the resolution, after continuum extrapolation, were 
high enough that the lattice and pcrturbativc curves could be subtracted from each other 
(rather than normalizing them to a function which diverges at short distances) , then the two 
features could be disentangled by inspecting intermediate distances, 0.2 < r < 0.3, in which 
the threshold region contributes much more prominently than the transport peak. In contrast 
the two features are difficult to tell apart if reliable data is only available for 0.3 < r < 0.5. 
(That said, we stress again that the present section is meant as indicative only.) 



7. Conclusions 

The main purpose of this paper has been to compute the massive quark-number susceptibility 
and vector current correlator at next-to- leading order (NLO) in thermal QCD. The NLO 
results are shown in eqs. (4.1), (4.4), (4.5), and illustrated numerically in figs. 1, 2. 

Our semi-analytic results can be directly compared with numerical lattice Monte Carlo 
simulations. Although no continuum limit has been reached for the massive vector current 
correlator [21], the agreement as seen in figs. 1, 2 is quite remarkable. (We have compared 
with quenched data, because these are "closest" to the continuum limit, but our analytic 
results are also valid for the unquenched case.) 

The good agreement suggests that resummed perturbative NLO computations, through 
which the Minkowskian spectral function has been determined around the threshold re- 
gion [51, 52], might be more accurate than sometimes assumed. A similar observation has also 

^The approach in this section is schematic; for possible other model shapes see, e.g., ref. [53]. (A general 
discussion of the sum rule approach can be found in ref. [54].) 
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been made by comparing NLO computations to lattice data at spatial separations relevant 
for quarkonium physics [55]. 

The back side of a remarkable agreement is that the strict NLO result of the present paper 
reflects the physics of an infinitely narrow transport peak, rather than genuine heavy quark 
diffusion, and of a threshold singularity, rather than genuine quarkonium resonances — and 
yet it works well. This implies that even though impressive efforts to extract properties of 
the transport peak and quarkonium resonances from lattice data are being made, a very fine 
resolution is needed for getting systematic errors fully under control. 

In order to quantify the resolution that lattice simulations should reach for observing 
substantial deviations from the NLO results, we have finally probed the significance of possible 
non-perturbative effects. The impact of a smeared transport peak is illustrated in fig. 3; the 
impact of resonance-like spectral weight below the threshold in fig. 4. Interpreting lattice 
results such as those in rcf. [21] as a deviation from NLO, it appears that a diffusion coefficient 
2ttTD ~ 1 — 2 could be observed in principle, but that larger values arc hard to disentangle. 
Perhaps, a fruitful approach is to determine the diffusion coefficient by other means [14, 15], 
subtract the corresponding contribution from lattice data, and hope that any remaining 
deviations reflect the physics of the threshold region. 

Ultimately, going beyond modelling, we suspect that systematic studies are only possible 
once a continuum extrapolation is available in a broad r-range. At this stage vacuum physics 
can be subtracted, as outlined in rcf. [27], and the non-divergent remainder subjected (at 
least in principle) to a model- independent analytic continuation [30]. 
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Appendix A. Results for individual master sum-integrals 

In this appendix we list results for the 2-loop sum- integrals in eqs. (2.13), (2.17), after carrying 
out the Matsubara sums as well as the Fourier transformation in eq. (2.18). For brevity we 
refer to the sum-integrals with the notation 



■j-mim2ms 



l^n K{P} 



K-Q 



(A.l) 



Q={u„,0) 
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The index mi guarantees that all the masters have the same dimensionality. It is also 
convenient to introduce the shorthand notations 



Aar = ek + (TEp + TEpk 



Ep + a Epk , 



(A.2) 



where the energies are defined as in cq. (2.11). In some cases the remaining integrands can 
be simplified via the symmetrization p -f-T- k — p. 

Four of the masters appearing are independent of the external momentum Q, and therefore 
lead to contact terms: 
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(A.6) 



P^P-K 



No further details are given since these contact terms cancel against contributions from 
masters with m2 > (eqs. (A.9), (A. 13), (A. 16), (A. 18), (A.24)). After carrying out the 
Matsubara sums the remaining masters can be represented as (with the notation of eq. (3.1)) 
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Appendix B. Renormalization of the complete result 

Inserting the expressions for the master sum-integrals from appendix A into eqs. (2.12), 
(2.13), (2.15), (2.17), we obtain renormahzed results for the correlators considered. Like in 
sec. 3, the result can be divided into r-dependent and constant terms. The former reads 



^NLO I 

V I r-dep. 



+ 



+ 



+ 



D 



(r) 



p,k ^kEpEpkA_^__A__^_ 
Jp,k ^kEpEpkA_^__A_^ 

f 2-^efc-Epfc 

Jp,k ^kEpEpkA_^__^_A — 
/ ^2E,(T)|"eq. (B.4)" 



-Ei 



-Ei 



E^k + M^ 



E^pk + M^ 



+El + E%-M^ 



A 



+ 



Ah-A_+, 



■++ 



A. 

'A, 



+ 



24 



A, A 



■+ 



A 



+ 



A++A__ 



+ 



+ 



A^+A^._A 



A2 A, 



-A-+ 



A?-+A++A__ 



19 



+ 



+ 



riBiek) 



1 m2 M2(2£'2 + M2)/ 1 



^2 + £;4 + 
p p 



pk 

^p^pk 



1 



+ 



£■4 



+ 



-+/ J 



+ 



npiEpk) 



P 



pk 



1 M2 M2(2£2^M2)/ 1 



2i?2 + i5;4 + 



'^{El+El^-M 



En 



1 



+ 



+ 



E'p [ 



A^A++A_ 
1 



A+A++A_ 



J^EpdE^D,EAT)\0 

nsiek) 



n^iEpk) 



El 2E^ 



2E^ 



\Al^ ' A2 
1 

A_A+_A_ 
1 

' A_A+_A_ 



+ 



1 



1 



A^ A^, 



2£;|A+A_ 



A++A+_ 



E, 



pk 



1 ^ M'^(2El + M'^) ^ l_ 



El 2^4 



A__A_ 
1 



2EI 



VA+H-A__ 



+ 



A+-A_, 



(B.l) 



The constant contribution, in turn, can be expressed as 
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Finally, for the susceptibility, the NLO term amounts to 
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The "0"s in eqs. (B.l), (B.2), (B.3) represent vacuum contributions that vanish after renor- 
malization. The coefficients of d^^D2£^{T) and Tn'l{Ep) vanish already when eqs. (2.12), 
(2.13) are summed together, but the coefficients of D2Ep{'T) and Tn'^{Ep) not. They can be 
expressed as (a principal value prescription is implied where necessary) 
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After having expressed the integrals in these forms, we can make use of Lorentz invariance 
in order to remove redundant structures: 
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Inserting these relations, eq. (B.5) vanishes exactly. The coefficient of D2Ep{t), eq. (B.4), 
does not vanish; after the transformation of eqs. (B.17)-(B.20) it can be written as 
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where the first row represents the result of the ultraviolet sensitive integrals containing single 
or double propagators. The remaining k-intcgral is infrared (IR) divergent and needs to be 
evaluated together with the other terms of eq. (B.l). Its form after angular integration can 
be found in eq. (4.4), as the structure preceding 9{k). 

It may be interesting to note that regulating the IR through a "gluon mass" A, we obtain^ 
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Here P* = {~iEp, p). The sum of these (with the proper coefficients) amounts to a particular 
choice of j^dk6{k)/k, appearing on the first row of eq. (4.4); a corresponding subtraction 
of 9{k) on the sixth row of eq. (4.4) would remove all powerlike terms at large k, with the 
price of inserting more structure into the infrared regime k ^ X. 
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